home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 476-500 / disk_499 / diglib / diglib.lzh / source / CONTOR.for < prev    next >
Text File  |  1991-04-13  |  7KB  |  223 lines
  1.         SUBROUTINE CONTOR(Z,NZ,IZ,MX,MY,X1,XMX,Y1,YMY,NL,CL)
  2.         IMPLICIT NONE
  3. C
  4. C       THIS SUBROUTINE WILL PRODUCE A CONTOUR PLOT OF THE FUNCTION
  5. C       DEFINED BY Z(I,J) = F(X(I),Y(J)).   IT IS ASSUMED THAT
  6. C       A CALL TO "MAPIT" HAS ALREADY BEEN MADE TO ESTABLISH THE
  7. C       COORDINATE AXIS (X,Y), WITH X LIMITS COVERING THE RANGE
  8. C       X1 TO XMX, AND Y LIMITS COVERING THE RANGE Y1 TO YMY.
  9. C
  10. C       Bad bug fixed: 23-APR-1985 by Hal R. Brand
  11. C
  12. CArguments:
  13. C
  14. C  Input
  15. C
  16. C       Z               * Type: real array.
  17. C                       * The values of the function to contour:
  18. C                          Z(I,J) = F(Xi,Yj) where:
  19. C                            Xi = X1 + (i-1)*(XMX-X1)/(MX-1)
  20. C                            Yj = Y1 + (j-1)*(YMX-Y1)/(MY-1)
  21. C
  22. C       NZ              * Type: integer constant or variable.
  23. C                       * The first dimension of the array Z - not necessaril
  24. C                               equal to MX, but MX <= NZ.
  25. C
  26. C       IZ              * Type: byte array.
  27. C                       * Used internally for working storage.
  28. C
  29. C       MX              * Type: integer constant or variable.
  30. C                       * The number of X grid points.
  31. C
  32. C       MY              * Type: integer constant or variable.
  33. C                       * The number of Y grid points.
  34. C
  35. C       X1              * Type: real constant or variable.
  36. C                       * The minimum X value.
  37. C
  38. C       XMX             * Type: real constant or variable.
  39. C                       * The maximum X value.
  40. C
  41. C       Y1              * Type: real constant or variable.
  42. C                       * The minimum Y value.
  43. C
  44. C       YMY             * Type: real constant or variable.
  45. C                       * The maximum Y value.
  46. C
  47. C       NL              * Type: integer constant or variable.
  48. C                       * The number of contour levels.
  49. C
  50. C       CL              * Type: real array.
  51. C                       * The coutour levels to draw.   (Same units as
  52. C                               F() or Z().)
  53. C
  54. C  Output
  55. C
  56. C    None.
  57. C
  58. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  59. C
  60.         REAL*4 Z(NZ,MY),X1,Y1,YMY
  61.         REAL*4 CL(NL)
  62.         INTEGER*1 IZ(MX,MY)
  63.         REAL*4 CLEVEL,XL,DX,YL,DY
  64.         INTEGER IOLD,JOLD,IN,JN,NX,NY,NLOOP
  65.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  66.      1   NX,NY,XL,DX,YL,DY
  67.         INTEGER I,J,NL,ICIR,IU,JU,NC
  68. C
  69. C       INITIALIZE ROUTINE
  70. C
  71.         XL = X1
  72.         YL = Y1
  73.         DX = XMX-X1
  74.         DY = YMY-Y1
  75.         NX=MX
  76.         NY=MY
  77.         NLOOP=MIN1(FLOAT(NX)/2.0+.5,FLOAT(NY)/2.0+.5)
  78. C       START SEARCHING FOR PLUS-MINUS TRANSITIONS
  79. C       TO START A CONTOR ON.
  80.         DO 50 NC=1,NL
  81. C
  82. C       ZERO ARRAY SHOWING WHERE WE HAVE BEEN
  83. C
  84.         DO 2 J=1,NY
  85.                 DO 1 I=1,NX
  86.                         IZ(I,J)=0
  87. 1                       CONTINUE
  88. 2               CONTINUE
  89.         CLEVEL=CL(NC)
  90.         DO 50 ICIR=1,NLOOP
  91.                 IU=NX+1-ICIR
  92.                 JU=NY+1-ICIR
  93.                 DO 10 J=ICIR,JU-1
  94.                         CALL LOOK(Z,ICIR,J,1,IZ,NZ,NX)
  95. 10                      CONTINUE
  96.                 DO 20 I=ICIR,IU-1
  97.                         CALL LOOK(Z,I,JU,2,IZ,NZ,NX)
  98. 20                      CONTINUE
  99.                 DO 30 J=JU,ICIR+1,-1
  100.                         CALL LOOK(Z,IU,J,3,IZ,NZ,NX)
  101. 30                      CONTINUE
  102.                 DO 40 I=IU,ICIR+1,-1
  103.                         CALL LOOK(Z,I,ICIR,4,IZ,NZ,NX)
  104. 40                      CONTINUE
  105. 50              CONTINUE
  106.         RETURN
  107.         END
  108. C
  109. C
  110. C
  111.         SUBROUTINE LOOK(Z,II,JJ,M,IZ,NZ,IZDIM)
  112.         INTEGER*1 IZ(IZDIM,2)
  113.         DIMENSION Z(NZ,2)
  114. C
  115. C       This subroutine looks for a contour starting at the point (II,JJ)
  116. C       with the contour being oriented such that the point (II,JJ) is
  117. C       greater than the current contouring level, and its neighbor (specifie
  118. C       by M) is less than the current contouring level.
  119. C
  120.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  121.      1   NX,NY,XL,DX,YL,DY
  122.         DIMENSION IDMODE(3,4)
  123.         DATA IDMODE/4,1,2,  1,2,3,  2,3,4,  3,4,1/
  124. C
  125.         IOLD=II
  126.         JOLD=JJ
  127.         MODE=M
  128.         CALL NEWP(1,MODE)
  129. C
  130. C       LOOK FOR CONTOR STARTING HERE.   THE "OLD" POINT IS ALWAYS THE
  131. C        POSITIVE ONE, SO THE TEST IS EASY.
  132. C
  133.         IF (Z(IOLD,JOLD) .GE. CLEVEL .AND. Z(IN,JN) .LT. CLEVEL) GOTO 20
  134. 100     RETURN
  135. C
  136. C       CHECK FOR CONTOR PREVIOUSLY THRU HERE - "SEGMNT" RETURNS THE POINT
  137. C        WE MARK WHEN EVER A CONTOUR PASSES THRU THE POINTS "(IOLD,JOLD)"
  138. C        AND "(IN,JN)".   "SEGMNT" ALSO RETURNS THE MARK THAT SHOULD BE
  139. C        PLACED GIVEN THE ORIENTATION OF THE CONTOUR.
  140. C
  141. 20      CALL SEGMNT(ICI,ICJ,ISEG)
  142.         IF ((IZ(ICI,ICJ).AND.ISEG) .NE. 0) RETURN
  143. C
  144. C       NEW CONTOUR.   TRACE IT TILL IT ENDS BY LOOPING BACK ON ITSELF, OR
  145. C        RUNNING OFF THE GRID.
  146. C
  147.         CALL ZPNT(XX,YY,Z,NZ)
  148.         CALL SCALE(XX,YY,VX,VY)
  149.         CALL GSMOVE(VX,VY)
  150.         IOLD=IN
  151.         JOLD=JN
  152. 30              CONTINUE
  153.                 DO 50 N=2,4
  154.                         CALL NEWP(N,MODE)
  155.                         IF (IN .LT. 1 .OR. IN .GT. NX) GO TO 100
  156.                         IF (JN .LT. 1 .OR. JN .GT. NY) GO TO 100
  157.                         IF (SIGN(1.0,Z(IOLD,JOLD)-CLEVEL) .NE.
  158.      1            SIGN(1.0,Z(IN,JN)-CLEVEL)) GO TO 60
  159.                         IOLD=IN
  160.                         JOLD=JN
  161. 50                      CONTINUE
  162. C
  163. C       IT IS IMPOSSIBLE TO JUST FALL THRU DUE TO THE ALGORITHM
  164. C
  165. C               STOP 'SERIOUS ERROR'
  166. 60              CONTINUE
  167. C
  168. C       FOUND THE NEXT INTERSECTION.   SEE IF IT HAS ALREADY BEEN MARKED.
  169. C        IF SO, THEN WE ARE DONE, ELSE, MARK IT, DRAW TO IT, AND CONTINUE ON.
  170. C
  171.                 CALL SEGMNT(ICI,ICJ,ISEG)
  172.                 IF ((IZ(ICI,ICJ).AND.ISEG) .NE. 0) RETURN
  173.                 IZ(ICI,ICJ)=(IZ(ICI,ICJ).OR.ISEG)
  174.                 CALL ZPNT(XX,YY,Z,NZ)
  175.                 CALL SCALE(XX,YY,VX,VY)
  176.                 CALL GSDRAW(VX,VY)
  177.                 MODE=IDMODE(N-1,MODE)
  178.                 GO TO 30
  179.         END
  180. C
  181. C
  182. C
  183.         SUBROUTINE SEGMNT(ICI,ICJ,ISEG)
  184.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  185.      1   NX,NY,XL,DX,YL,DY
  186.         ICI=MIN0(IOLD,IN)
  187.         ICJ=MIN0(JOLD,JN)
  188.         ISEG=1
  189.         IF (IOLD .EQ. IN) ISEG=2
  190.         RETURN
  191.         END
  192. C
  193. C
  194. C
  195.         SUBROUTINE NEWP(I,M)
  196.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  197.      1   NX,NY,XL,DX,YL,DY
  198.         DIMENSION IDELI(4),JDELJ(4)
  199.         DATA IDELI,JDELJ / 0,1,0,-1,   1,0,-1,0/
  200.         INDEX=MOD(2+I+M,4)+1
  201.         IN=IOLD+IDELI(INDEX)
  202.         JN=JOLD+JDELJ(INDEX)
  203.         RETURN
  204.         END
  205. C
  206. C
  207. C
  208.         SUBROUTINE ZPNT(X,Y,Z,NZ)
  209.         DIMENSION Z(NZ,2)
  210.         COMMON /CONTR/ CLEVEL,IOLD,JOLD,IN,JN,
  211.      1   NX,NY,XL,DX,YL,DY
  212.         A=Z(IN,JN)-Z(IOLD,JOLD)
  213. C       IF NO CHANGE IN Z'S, PICK OLD POINT SO AS TO STAY TO RIGHT
  214.         IF (A .EQ. 0.0) GO TO 10
  215.         A=(CLEVEL-Z(IOLD,JOLD))/A
  216. 10      X=A*(IN-IOLD)+IOLD
  217.         Y=A*(JN-JOLD)+JOLD
  218. C       NOW CONVERT INDEXS TO X,Y VALUES
  219.         X=(X-1.0)*DX/(NX-1)+XL
  220.         Y=(Y-1.0)*DY/(NY-1)+YL
  221.         RETURN
  222.         END
  223.